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ABSTRACT 

Context. When massive stars exert a radiation pressure onto their environment that is higher than their gravitational attraction (super- 
Eddington condition), they launch a radiation-pressure-driven outflow, which creates cleared cavities. These cavities should prevent 
any further accretion onto the star from the direction of the bubble, although it has been claimed that a radiative Rayleigh-Taylor 
instability should lead to the collapse of the outflow cavity and foster the growth of massive stars. 

Aims. We investigate the stability of idealized radiation-pressure-dominated cavities, focusing on its dependence on the radiation 
transport approach used in numerical simulations for the stellar radiation feedback. 

Methods. We compare two different methods for stellar radiation feedback: gray flux-limited diffusion (FLD) and ray-tracing (RT). 
Both methods are implemented in our self-gravity radiation hydrodynamics simulations for various initial density structures of the 
collapsing clouds, eventually forming massive stars. We also derive simple analytical models to support our findings. 
Results. Both methods lead to the launch of a radiation-pressure-dominated outflow cavity. However, only the FLD cases lead 
to prominent instability in the cavity shell. The RT cases do not show such instability; once the outflow has started, it precedes 
continuously. The FLD cases display extended epochs of marginal Eddington equilibrium in the cavity shell, making them prone 
to the radiative Rayleigh-Taylor instability. In the RT cases, the radiation pressure exceeds gravity by 1-2 orders of magnitude. The 
radiative Rayleigh-Taylor instability is then consequently suppressed. It is a fundamental property of the gray FLD method to neglect 
the stellar radiation temperature at the location of absorption and thus to underestimate the opacity at the location of the cavity shell. 
Conclusions. Treating the stellar irradiation in the gray FLD approximation underestimates the radiative forces acting on the cavity 
shell. This can lead artificially to situations that are affected by the radiative Rayleigh-Taylor instability. The proper treatment of direct 
stellar irradiation by massive stars is crucial for the stability of radiation-pressure-dominated cavities. 
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1. Introduction 

Outflows are a unique characteristic of the star formation pro- 
cess in general. The strong radiation pressure provided by mas- 
sive luminous stars plays an important role in the dynamics 
of these outflows. Observations show that massive star-forming 
regions involve several large-scale outflows and their interac- 
tion with the interstellar material depicts an important part of 
the complex morphology of the cluster center (e.g. Shepherd 
& Churchwell 1996; Zhang et al. 2001; Beltran et al. 2006; 
Zhang et al. 2007; Fallscheer et al. 2009; Beuther et al. 2010; 
Quanz et al. 2010). Non-radiative MHD simulations (Banerjee 
& Pudritz 2006, 2007; Hennebelle el al. 201 1) have shown that 
cavities of lower density potentially form before the onset of star 
formation during the collapse of magnetized clouds. The treat- 
ment of the radiation field is the next step in investigating - es- 
pecially the thermodynamics of - the cavities. 

From the perspective of numerical development, the imple- 
mentation and improvement of radiation transport schemes in 
(magneto-)hydrodynamics simulations is both a challenging and 
timely problem in modern astrophysics, especially in the field of 
(massive) star formation, in which the radiative feedback onto 
the environment plays a crucial role (Yorke & Sonnhalter 2002; 
Krumholz et al. 2007, 2009; Price & Bate 2009; Commercon 
et al. 2010; Bate 2010; Kuiper et al. 2010a, 201 1). 

The hybrid radiation transport scheme applied in this pa- 
per was previously used in our studies of the radiation pres- 



sure barrier in the formation of massive stars (Kuiper et al. 
2010a) and the angular momentum transport in disks around 
massive stars (Kuiper et al. 201 1). We introduced the derivation 
of the solving algorithm as well as its numerical implementation 
into a three-dimensional (3D) (magneto-)hydrodynamics code 
in Ktiiper et al. (2010b). The solver algorithm superposes two 
radiation fields: a frequency-dependent ray-tracing (RT) compo- 
nent representing the stellar irradiation field and a frequency- 
averaged (gray) flux-limited diffusion (FLD) component repre- 
senting the thermal dust emission. These splitting schemes were 
previously used in one-dimensional (ID) simulations ( vVoliire 
& Cassinelli 1986, 1987; Edgar & Clarke 200 ) and in two- 
dimensional (2D) hydrostatic disk atmosphere models (Murray 
et al. 1994) . 

The RT methods are commonly used in radiative trans- 
fer models without hydrodynamical motion (e.g. Efstathiou 
& Rowan-Robinson 1990; Steinacker et al. 2003). Alternative 
methods for these kinds of radiative transfer models include 
the short characteristics method (DuUemond & Turolla 2000) 
and Monte Carlo radiative transfer (Bjorkman & Wood 2001; 
DuUemond & Dominik 200-I-). Comparison benchmark studies 
of these methods were presented in Pascucci et al. (2004) and 
Pinte et al. (2009). These methods are able to solve the radia- 
tion transport problem to high accuracy for a fixed static con- 
figuration, but they require much CPU time, which yields low 
efficiency in time-dependent hydrodynamics simulations. 
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In this sense, the FLD approximation (Kley 1989; 
Bodenheimer et al. 1990; Kiahr et al. 1999; Klahr & Kley 2006) 
provides a fast method to determine the temperature evolution 
accurately in the optically thick (diffusion) and optically thin 
(free-floating) limit. However, the approximation becomes in- 
accurate in multi-dimensional anisotropic problems as well as 
in transition regions from optically thin to thick (shown e.g. in 
Boley et al. 2007). One example of this transition region is the 
cavity shell investigated in this paper The FLD approximation 
in the frequency-averaged (gray) limit remains the default tech- 
nique in modem multi-dimensional radiation hydrodynamics. 
In an exceptional case, a frequency-dependent FLD solver was 
used in Yorke & Sonnhalter (2002). 

From the theoretical point of view, the stability of radiation- 
pressure-dominated cavities has a direct impact on the star for- 
mation efficiency and accretion onto a luminous massive star 
Krumholz et al. (2009) proposed that massive stars can grow 
beyond their Eddington limit because these radiation-pressure- 
dominated cavity shells are subject to the so-called "radia- 
tive Rayleigh-Taylor instability" allowing additional mass to be 
fed onto the central accretion disk. Following the approach of 
Nakano (1989) and Yorke & Sonnhalter (2002), we proposed 
the feeding of a massive star beyond its Eddington limit by disk 
accretion only (Kuiper et al. 2010a, 2011). While the anisotropy 
of the thermal radiation field sufficiently diminishes the radiation 
pressure onto the disk accretion flow near the midplane, the radi- 
ation pressure in the polar direction is able to launch an outflow, 
forming a stable cavity. 

The removal of a substantial fraction of the initial core mass 
by the radiative launching of outflows influences the star for- 
mation efficiency in the central core region. These feedback ef- 
fects are proposed to play an important role in explaining the low 
star-formation efficiencies observed (McKee & Ostriker 2007). 
Simulations of Hll regions in massive star forming regions by 
Dale & Bonnell (2011) and close to massive stars forming in 
the cluster center by Peters et al (2010) suggest that ionized 
gas fills preexisting voids and bubbles. Long-lived cleared polar 
cavities therefore would make the expansion of HII regions and 
their interaction with the stellar cluster environment easier. The 
feedback of the energy injection from jets and outflows onto the 
stellar cluster formation is briefly discussed in Bonnell & Bate 
(2006). 

In this study, we investigate different implementation meth- 
ods for direct stellar irradiation feedback (see Sect. 3) and their 
influence on the stability of radiation pressure dominated cavi- 
ties. We present in Sect. 4 the qualitative outcome and Sects. 5 
and 6 the quantitative analyses of the simulations. The observed 
difference in the radiative acceleration of the cavity shell depend- 
ing on the applied radiation transport method is analytically de- 
rived in Sect. 7. Finally, a comprehensive comparison section 
(Sect. 8) to previous simulations, analytic work, and observa- 
tions as well as a brief summary (Sect. 9) are provided. 

2. Physical and numerical model 

2.1. Physics included 

We compare the effect of two different radiation transport meth- 
ods on stellar radiation feedback. The goal is to investigate 
the stability of the shell surrounding the radiation-pressure- 
dominated cavity. We emphasize that in these simulations we do 
not aim to provide a complete description of an outflow region 
around a massive star. 

A list of potential limitations and caveats includes: 



1. The coUimation effects by magnetic fields - which are cer- 
tainly important for the morphology of the cavity - are ne- 
glected. Recent results of MHD simulations of jet formation 
are summarized by Fendt et al. (201 1). Potentially, the rela- 
tive importance of radiative or magnetic forces in jet forma- 
tion can be distinguished observationally based on morphol- 
ogy: a magnetically launched fast jet would be more colli- 
mated than a radiation-pressure-driven wide-angle outflow 
(see e.g. Arce et al. 2007; Pudritz et al. 2007, and citations 
therein). 

2. Another example of the influence of magnetic fields, the so- 
called photon bubble instability (rurner el al. 2007) could 
potentially diminish the radiation pressure and could there- 
fore alter the morphology of a cavity filled with magnetized 
gas. 

3. Evaporation of dust is included, but otherwise the dust is as- 
sumed to be perfectly coupled to the gas. 

4. The dust opacity clearly dominates the absorption in any 
dusty region, but in the dust-free zone around the massive 
star up to the dust sublimation front the (remaining) gas 
opacity is set to be constant /Cgas = 0.01 cm^ g"' . Such a con- 
stant gas opacity of course denotes only a first order approx- 
imation to the complex molecular line absorption features. 
An optically thick gas disk could e.g. amplify the flashlight 
effect ( ianaka & iNakamolo 201 1) and therefore increase the 
radiative flux into the polar regions. 

5. Since the cavity growth is driven by radiation pressure, the 
quantitative details of the launching process and the growth 
phase depend on the stellar evolution model as well as the 
dust model. Future studies will attempt to establish the de- 
tails of this dependence. 

6. The low-density cavity is potentially enriched, but certainly 
influenced, by a stellar wind from the central massive star, 
which is neglected herein. Early proto-stellar winds could, 
e.g., lead to an early formation of rather collimated polar 
cavities ( ^unningham et al, 201 1). The radiation pressure at 
later epochs in such a configuration could be channeled into 
pre-existing cavities, leading to smaller opening angles for 
the radiation pressure dominated cavities. 

7. In Peters el al. (2010, 201 ), ionizing radiation was proposed 
to dominate the cavity dynamics in regions of massive star 
formation. 

2.2. Numerical resolution 

To obtain unambiguous results, the relevant physical length 
scales have to be resolved by the numerical simulations. 
Therefore, we discuss these scales with respect to the resolution 
of our radiation hydrodynamical simulations. We demonstrate 
that the simulations resolve all relevant length scales sufficiently. 

A study of a radiative Rayleigh-Taylor instability in the shell 
surrounding the outflow cavity needs to resolve the wavelengths 
of this instability along the shell discontinuity. At small wave- 
lengths, the radiative Rayleigh-Taylor instability is likely to be 
suppressed by diffusion, but at large wavelengths - compara- 
ble to the physical size of the cavity - the instability occurs. 
Quantitatively, the shortest unstable wavelength for instance is 
determined to be 1000 AU for a IOMq and 10000 AU for a 
100 Mq star ( , , Krumholz 201 1). The polar resolution 
rA6 of the spherical grid in our simulations grows linearly with 
the radius and is fixed to rA6 < O.lr. The grid size along the shell 
discontinuity is typically a factor of ten smaller than the short- 
est unstable wavelength. Hence, the length scale of the radiative 
Rayleigh-Taylor instability is fully resolved. 
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Fig. 1. The cooling length scale (dashed line) of the thermal 
dust emission in the gray approximation as a function of the ac- 
tual location of the cavity shell /^cavity The solid line denotes the 
resolution of our numerical grid. 



Two other important length scales are related to the radiative 
properties of the shell. The "cooling length scale" k is given by 
an optical depth of t = 1 in the long wavelength regime, and 
the "irradiation length scale" h is given by an optical depth of 
T = 1 of the broad stellar irradiation spectrum. Fig. 1 shows the 
cooling length scale in the gray approximation as a function of 
the location of the cavity shell /^cavity The shell density is chosen 
to be p = 4 X 10"'^ g cm"^ as depicted in Fig. 5. The opacities 
are computed as the Rosseland mean opacities of Laor & Draine 
(1993). The stellar temperature is taken from the Hosokawa & 
Omukai (2009) tracks for a 20 Mq star and an accretion rate of 
M - 10"^ Mg yr"', and the temperature at the location ^cavity 
is computed with a slope T cc R^^^^^^ in the gray approximation. 
For simplification, the evaporation of dust grains is neglected, 
but would result in an even larger cooling length-scale. Hence, 
we can be sure that the length scale of cooling is resolved. 

The absorption length scale Z» of the stellar irradiation spec- 
trum is much smaller than the cooling length scale for ther- 
mal dust emission. Using again the opacities of Laor & Draine 
(19* ), Fig. 2 shows the absorption length scale of the whole 
stellar irradiation spectrum for different gas densities. For the 
aforementioned shell density ofp - 4x 10"'^ g cm"-* and a typ- 
ical stellar temperature at the interesting point in time at the po- 
tential onset of instability of about 10000 K < < 30000 K, the 
absorption length scale of the photons with frequencies higher 
than the peak of the stellar black-body spectrum is about 10 to 
100 AU; these photons are most likely to be absorbed in the 
first grid cell on top of the cleared cavity. For lower frequency 
photons, the absorption length scale smoothly rises up to the 
1000 AU scale and is most likely to be resolved. In the long 
wavelength regime of the spectrum, the results of the previous 
paragraph for the thermal dust emission of course apply. 

The morphology in the radial direction during the cavity 
growth phase can be distinguished into three important regimes: 
the inner low-density cleared cavity, the shell of swept-up mate- 
rial on top of this cavity, and the remnant of the in-falling enve- 
lope on larger scales. The resolved transition of the gas density, 
temperature, and velocity in these regimes is depicted in the con- 
text of the following analyses in Figs. 3 (two-dimensional snap- 
shots) and 5 (along the polar axis). Even the steep gradient of 
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Fig. 2. The absorption length scale h regarding the stellar irra- 
diation spectrum as a function of frequency v or wavelength A, 
respectively. From top to bottom, the different solid lines denote 
a density of p = lO^''^ IQ-^^, 1Q-" , 10-'^ and lO^'^ g cm-\ 
respectively. The two horizontal lines mark the highest and 
lowest resolution of the spherical grid. The right vertical axis 
shows the scale for the stellar spectra: From top to bot- 
tom, the dashed lines denote black body spectra of r» = 
100000,30000, 10000, and 3000 K. The effect of potential dust 
evaporation is not included here, but would lead to a higher reso- 
lution of the numerical grid in terms of the UTadiation spectrum. 



the density at the inner rim of the cavity shell (Fig. 5, upper left 
panel) is resolved by four grid cells, and the smooth decrease to 
larger radii by dozens of grid cells. Therefore, the shell morphol- 
ogy is clearly resolved in the simulations. 

2.3. Axial symmetry 

The simulations in this paper are performed assuming axial sym- 
metry. We are convinced that this is not a critical limitation of 
the result we present. In 2D simulations with the same initial 
conditions, but different methods for stellar radiation feedback, 
we observe strong differences in the resulting radiation feedback 
onto the shell on top of the cleared cavity. These differences do 
not depend on the detailed morphology of this shell and there- 
fore the conclusion is also applicable to non-axially symmetric, 
3D configurations. 

In our 2D simulation results with the FLD approximation, 
the cavity shell undergoes an instability; this result matches the 
3D simulation result of Krumholz et al. (2009). In our 2D simu- 
lations as well as in our 3D simulation (Kuiper et al. 201 1 ), in- 
cluding the RT step for direct irradiation feedback, the radiation- 
pressure-dominated cavity remains stable. 

3. Method 

3.1. Equations 

To follow the motion of the gas, we solve the equations of com- 
pressible self-gravity hydrodynamics, including shear-viscosity 
as described in Kuiper et al. (2010a). For this task, we use the 
open source MHD code Pluto (Mignone et al. 2007). This set 
of equations is coupled to the radiation transport in the medium. 
In contrast to our previous studies, we investigate the influence 
of two different methods to determine the radiative flux, namely 
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our hybrid radiation transport scheme (ray-tracing + flux-limited 
diffusion) and flux-limited diffusion only. 

3.1.1. Ray-tracing -i- flux-limited diffusion 

For simulations labeled "RTh-FLD", we use our hybrid radia- 
tion transport module (see Kuiper el al. 20101 ). The total ra- 
diation energy density Ftot is divided into the ffux f,(v, r) of 
the frequency-dependent stellar irradiation and the ffux of the 
frequency-averaged thermal radiation energy density F 



d,E^+f,V-F = -/,(V-f.-e+), 



(1) 
(2) 



where Eq. (1) denotes the evolution of the thermal radiation en- 
ergy density /Jr. The factor /c = (cyp/A-aT^ + l) depends 
only on the ratio of internal to radiation energy and contains 
the specific heat capacity, Cy, and the radiation constant, a. The 
source term depends on the physics included and contains 
any additional energy source such as hydrodynamical compres- 
sion -PV ■ u and viscous heating. We solve Eq. (1) using the 
so-called ffux-limited diffusion approximation, in which the ffux 
is set proportional to the gradient of the radiation energy den- 
sity (F = -ZJVEr). The diffusion constant D = AclpK^ depends 
on the ffux limiter A and the Rosseland mean opacity /cr. The 
quantity c denotes the speed of light in a vacuum. We use the 
ffux limiter proposed by Levermore & Pomraning (1981) and 
neglect scattering. 

Eq. (2) calculates the ffux of the frequency-dependent stel- 
lar irradiation in a ray-tracing step. The first factor on the right 
hand side describes the geometrical decrease in the ffux propor- 
tional to r^^. The second factor describes the absorption of the 
stellar light as a function of the optical depth t(v, r) - K{v)p{r)r 
depending on the frequency-dependent mass absorption coeffi- 
cients a:(v). For this purpose, we use tabulated dust opacities by 
Laor & Draine (1993), including 79 frequency bins, and calcu- 
late the local evaporation temperature of the dust grains by using 
the formula of Isella & Nana (2005). The ffux at the inner radial 
boundary is given by the luminosity L», temperature r», and ra- 
dius R„ of the forming star For this purpose, we use tabulated 
stellar evolutionary tracks for accreting high-mass stars, calcu- 
lated by Hosokawa & Omukai (2009). The gas and dust temper- 
ature T is finally calculated in equilibrium with the combined 
stellar irradiation and thermal radiation field 



ar = ^R + 



kAT) c 



(3) 



with the Planck mean opacities k^. 

Numerical details, test cases, including a comparison of gray 
and frequency-dependent irradiation, as well as performance 
studies of the hybrid radiation transport scheme are summarized 
in Kuiper et al. (2010b). The viscosity prescription as well as 
the tabulated dust and stellar evolution model are presented in 
Kuiper et al. (2010a). 

3.1.2. Flux-limited diffusion only 

For simulations labeled "FLD", we do not ray-trace the stellar 
irradiation, but add the luminosity of the forming star as a source 
term to the right hand side of the FLD equation 



5,£R+/eV-f = /e (L.5(r) + e+) . 



(4) 



implemented as a Dirichlet boundary condition at the inner rim 
of the computational domain. The radiation energy at this inner 
boundary is therefore given by the stellar surface temperature 
and the assumption that the medium between the stellar surface 

and the inner rim at = 10 AU is optically thin. 

The hybrid radiation transport scheme (Sect. 3.1.1) as well 
as the FLD approach were tested in detail in comparison with 
a Monte Carlo solution using the RADMC code (Diillemond 
& Dominik 2004) or a short characteristics method using the 
RADICAL code (Dullemond & Turolla 20()i ') in a setup includ- 
ing a central star, an accretion disk, and an envelope as described 
in Pascucci et al. (2004) and Pinte et al. (2009). 

3.2. Numerical configuration 

The simulations are performed on a time-independent grid in 
spherical coordinates with a logarithmically stretched radial co- 
ordinate axis. The outer core radius is fixed to r^ax = 0.1 pc, 
and the inner core radius is fixed to rmin = 10 AU. The accu- 
rate size of this inner sink cell was determined in a parameter 
scan presented in Kuiper et al. (2010a), Sect. 5.1. The polar an- 
gle extends from 0° to 90° assuming midplane symmetry. The 
grid consists of 64 x 16 grid cells, i.e. the highest resolution of 
the non-uniform grid is chosen to be 



Ar X rA6 = 1.27 AU x 1.04 AU 



(5) 



Since the forming star is enclosed in the inner sink cell in our 
grid in spherical coordinates, the luminosity of the central star is 



around the forming massive star. The resolution decreases loga- 
rithmically in the radial outward direction proportional to the ra- 
dius. The radially inner and outer boundary of the computational 
domain are semi-permeable walls, i.e. the gas is allowed to leave 
the computational domain (by accretion onto the central star or 
due to radiative and centrifugal forces over the outer boundary), 
but cannot enter this domain. This outer boundary condition al- 
lows the mass reservoir for stellar accretion to be controlled by 
the initial choice of the mass of the pre-stellar core. 

The Pluto code uses high-order Godunov solver methods 
to compute the hydrodynamics, i.e. it uses a shock capturing 
Riemann solver within a conservative finite volume scheme. 
Our default configuration consists of a Harten-Lax-Van Leer 
Riemann solver and a so-called "minmod" ffux limiter us- 
ing piecewise linear interpolation and a Runge-Kutta 2 time- 
integration, also known as the predictor-corrector-method; for 
comparison we also refer to van Leer (1979). Therefore, the total 
difference scheme is accurate to second order in time and space. 

The internal iterations of the implicit solver of the FLD in 
Eqs. (1) or (4) is stopped at an accuracy of the resulting temper- 
ature distribution of either AT jT < 10"^ or Ar < 0.1 K. The 
internal iterations of the implicit solver for Poisson's equation is 
stopped at an accuracy of the resulting gravitational potential of 
A<I>/<D < 10 ^ 

3.3. initial conditions 

We start from a cold (Tq - 20 K) pre-stellar core of gas and 
dust. The initially constant dust-to-gas mass ratio is chosen to 
be Mdust/^gas = 0.01. Dynamically, the model describes a so- 
called quiescent collapse scenario without initial turbulent mo- 
tion (Mr - Ug -0) and the core is initially in slow solid-body ro- 
tation (\u^\/R = 0() = 5 * 10"'^ Hz). The total mass Mcore in the 
computational domain is chosen to be 50 or 100 Mq. In addi- 
tion to the inffuence of the different radiation transport methods 
described in the last sections, we check the dependence of the 
stability of the radiation pressure dominated cavity on the initial 
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density slope of the pre-stellar core. We chose the initial density 
slope to be either p oc r"'^ or p oc r"^. An overview of the runs 
is given in Table 1 . For a more comprehensive parameter scan of 



Run M,„,e pocr^ RTM Epoch of 



A-RT+FLD 


100 


-2 


RT-hFLD 


No 


A-FLD 


100 


-2 


FLD 


Yes 


B-RT-t-FLD 


50 


-2 


RT-l-FLD 


No 


B-FLD 


50 


-2 


FLD 


Yes 


C-RT+FLD 


100 


-1.5 


RT+FLD 


No 


C-FLD 


100 


-1.5 


FLD 


Yes 



Table 1. Overview of simulations performed. The table contains 
the run label (column 1), the two distinguishing parameters of 
the different initial conditions, namely the initial pre-stellar core 
mass Mcore (column 2) and the exponent y6 of the initial den- 
sity slope p(r) (column 3), in addition to the radiation transport 
method (RTM) in use (column 4), and the resulting existence of 
an epoch of marginal Eddington equilibrium E ^ I (column 5). 



the initial conditions of this pre-stellar core collapse model, we 
refer to Kuiper et al. (2010a) and Kuiper et al. (201 1). 

4. Qualitative simulation results 

We summarize the qualitative outcome, the stability, and the 
overall morphology of the radiation pressure dominated cavities, 
depending on the different initial conditions and the radiation 
transfer method. In the six simulations, the radiation pressure is 
generally high enough to launch an outflow and form a cleared 
polar cavity. 

In the case of the RT-i-FLD radiation transport method, these 
outflow cavities, once launched, rapidly grow in their extent up 
to the outer computational domain at 0. 1 pc away from the cen- 
tral star. In contrast to this monotonic growth, the outflow cav- 
ities in the FLD radiation transport scheme stop their increase 
along the polar axis at an extent of the order of several 100 to 
1400 AU depending on the initial conditions. This stopping of 
their growth is followed by a penetration of the gas mass on 
top of the outflow cavity formed so far, i.e. the outflow cav- 
ity collapses. During a follow-up second and third epoch of 
outflow launching, the resulting cavities grow in their extent, 
but on much longer timescales than in the RT-hFLD case. The 
epochs of frozen cavity growth in the FLD-only simulations 
suggest that (parts of) the top layer of the outflow cavity are 
in marginal Eddington equilibrium, i.e. the radiation pressure 
force is in equilibrium with the stellar gravity. The qualitative 
outcome of the simulations (depending on whether an extended 
epoch of marginal Eddington equilibrium occurs) is summarized 
in Table 1, column 5. 

In addition to the different growth rates, the radiation- 
pressure-driven outflow cavities also develop pronounced dif- 
ferences in their morphology depending on the applied radia- 
tion transport method. As an example, the different morpholo- 
gies produced by the different radiation transport methods are 
visualized in three snapshots in time during the onset of the out- 
flow launching in Fig. 3. In the case of FLD, the outflow is easily 
stopped by the infalling matter along the polar axis, while in the 
RT+FLD case the outflow cavity grows rapidly with time. In the 
FLD approximation, the radiative flux tends to follow a path that 
minimizes the optical depth and hence avoids the swept-up mass 



on top of the cavity. This avoidance is alleviated by the centrifu- 
gal forces, which diminish the gravitational attraction in regions 
slightly above the disk. In the ray-tracing method, the isotropic 
stellar irradiation flux directly impinges onto the swept-up mass 
on top of the polar cavity and pushes the mass to larger radii. 
The resulting large-scale morphology of the cavity is far more 
isotropic. The opening angle of the cavity is determined by the 
inner disk structure. 

In the following Sects. 5 and 6, we analyze quantitatively 
the outcome of the simulations presented here. We check our 
hypothesis of epochs of marginal Eddington equilibrium in the 
FLD-only runs and determine via analytical estimates (Sect. 7), 
why the cavity shells in the RT-hFLD case do not undergo these 
epochs and therefore remain stable with respect to the radiative 
Rayleigh-Taylor instability. 



5. Quantitative analysis of the cavity growth 

We analyze quantitatively the time-dependent extent of the out- 
flow cavity. The radial extent /^cavity of the outflow cavities above 
the central star is determined as the extent of the cleared cavity 
along the polar axis, as visualized in Fig. 3. The resulting cavity 
radii as a function of time are shown for the three different initial 
conditions as well as the different radiation transport methods in 
Fig. 4. 

For all initial conditions, the extent of the radiation pres- 
sure dominated cavities increases far more rapidly, if the stellar 
source of radiation is treated via a ray-tracing scheme (labeled 
"RT+FLD") than simply included in the flux-limited diffusion 
solver (labeled "FLD"). If using the FLD approximation, the 
outflow is launched a little bit earlier in time and after its ini- 
tial growth phase undergoes an epoch of marginal stability, in 
which the radiation pressure force along the polar axis seems 
to be balanced by gravity; the cavity growth along the polar 
axis is first stopped and then reversed. Under these conditions 
of marginal Eddington equilibrium (/rad ^ /grav) with FLD, the 
cavity shell region has been proposed to be subject to the ra- 
diative Rayleigh-Taylor instability ( .cquel & Krumholz 201 1). 
As a consequence, subsequent growth phases of the outflows in 
the FLD runs depend on the details of the subgrid models; the 
dynamics of the marginally balanced cavity shells are strongly 
influenced by the stellar evolution model and the dust model, 
particularly the treatment of sublimation and evaporation. 

In contrast to this epoch of marginal Eddington equilibrium, 
the extent of the outflow cavity in the RT + FLD simulations in- 
creases rapidly in time. With the exception of the Mcore = 50 M© 
case, in which a correspondingly lower mass star (yielding much 
lower luminosity) is formed, the outflow cavity increases in size 
monotonically in the RT + FLD simulations. 

The difference in the growth rate for the various radiation 
transport schemes is most prominent in the Mcore = 100 Mq and 
p oc case (Fig. 4 left panel), i.e. in the case of the initially 
highest mass in the central core region, allowing for a rapid for- 
mation and growth of the central star. In the Mcoie - 50 Mq case 
(Fig. 4 middle panel), the much lower stellar luminosity results 
in an initially unstable cavity region even in the RT+FLD case 
(but with a smaller extent than 100 AU). In the p oc r"'^ case 
(Fig. 4 right panel), the mass is initially more concentrated at 
larger radii and hence the cavity shell has to sweep up far more 
mass during its increase. Therefore, the growth phase of the cav- 
ity takes much longer than in the cases of initially steeper density 
profiles. 
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6. Quantitative analysis of tlie cavity sliell 
morpliology and dynamics 

To unveil the physical background of the qualitative and quan- 
titative difference of the outflow cavity structure in the simu- 
lations using either the RT-)-FLD or the FLD method, we now 
investigate the morphology and dynamics of the cavity shell de- 
pending on the radiation transport method. In this section, we 
analyze quantitatively the difference in the morphology and the 
dynamics of the cavity shell in the RT-hFLD and FLD cases. 
We focus on the data of the simulation with an initial core mass 
^coie = 100 Mq and a radial density slope of /3 = -2. 

The cavity in the FLD case increases in its first growth phase 
up to roughly 1400 AU before the expansion along the polar axis 
stops and the mass flow is reversed by gravity. We analyze the 
gas dynamics at the point in time when in both simulations (with 
RT+FLD and with FLD only) the expansion of the outflow cav- 
ity has arrived at the same location, namely at f = 16.7 kyr. The 
gas density, temperature, radial velocity, and the radial mass loss 
rate along the polar axis at this time are shown in Fig. 5. The 
upper left panel of this figure shows that there is almost no dif- 
ference in the radial extent, compactness, and peak density of the 
swept-up mass in the cavity shell. The peaks in the density distri- 
butions correspond to 4.1 x 10"'^ g cm"-' and 3.7 x 10"'^ g cm"-' 
in the RT-hFLD and FLD cases, respectively. 

The upper right panel of Fig. 5 shows that the RT-i-FLD 
run results in a continuously slightly higher temperature in the 
cleared cavity, the cavity shell, and beyond. The temperature dis- 
tribution in the pure FLD run stays 24% to 45% below the tem- 
perature distribution of the RT-i-FLD run. 

In contrast to these similarities, the two lower panels of Fig. 5 
highlight the striking difference of both runs in their radial ve- 
locities and mass flux along the polar axis. In the RT+FLD run, 
the swept-up shell material moves into the interstellar medium 
with a speed slightly higher than 100 km s"', which is roughly 
a factor of three higher than in the corresponding FLD run. At 
larger radii (r > 3000 AU), the velocity slope is still dominated 
by gravitational infall and hence is independent of the radiation 
transport method. 

Up to the onset of the Eddington equilibrium epoch in 
the FLD run, the mass of the central star in both simulations 
(RT-hFLD and FLD only) closely match, leading to the same 
strength of gravitational attraction for the matter in the cavity 
and shell towards the direction of the star Furthermore, the den- 
sity structure is roughly the same (Fig. 5, upper left panel) and 
the deviation in temperature within the cavity is smaller than 24- 
45% (Fig. 5, upper right panel), hence the slightly higher thermal 
pressure cannot be the main driver of the enormous difference in 
the radial velocities and mass flux (Fig. 5, lower panels). 

As a consequence, the radiation pressure acting on the swept- 
up mass on top of the outflow cavity has to differ in the two radi- 
ation transport methods. First, the FLD approximation assumes 
that the dust grains of the swept-up mass on top of the cavity 
are embedded in a local radiation field, whereas the RTh-FLD 
method takes into account the (frequency-dependent) absorption 
probability caused by direct stellar irradiation. The resulting dif- 
ference in the radiative acceleration at the top of the outflow cav- 
ity is analytically estimated in the following section. 

Secondly, since the FLD approximation is mathematically a 
moment method, the derivation of the FLD approximation in- 
cludes the integral over the angular distribution of the radiative 
flux; hence the emitted photons of the central star do not move 
along straight rays until they are absorbed. The radiative flux in 
the FLD approximation instead follows a path that minimizes 



the optical depth. The FLD method was originally introduced 
to describe spherically symmetric flows, hence the integral over 
the propagation direction did not result in unphysical behavior 
(see e.g. Bruenn 1985). In a multi-dimensional environment with 
a non-isotropic optical depth, this assumption of the FLD ap- 
proximation breaks down. Including higher-order terms in the 
derivation would minimize this inaccuracy and potentially yield 
a sufficient tracing of the correct photon path. However, in the 
RT+FLD scheme the irradiation by the central star is computed 
via a tay-tracing equation, which takes into account the correct 
photon propagation direction, i.e. the stellar irradiation flux di- 
rectly impinges on the swept-up mass shell at the top of the op- 
tically thin cavity. 

7. Analytical determination of the acceleration 
caused by stellar irradiation in the ray-tracing 
approach and the flux-limited diffusion 
approximation 

Our aforementioned analysis indicates that the difference be- 
tween the cavity shell stability in the two radiation transport 
methods is due to the higher radiative acceleration of the swept- 
up material when the RT approach is used for the stellar irradi- 
ation instead of the FLD approximation. The use of RT guaran- 
tees that the stellar radiation directly interacts with the swept-up 
mass at the top of the optically thin outflow cavity. The dust 
grains in the swept-up material therefore absorb the stellar flux 
following the frequency-dependent absorption coefficient a-(v). 
If we integrate over frequencies, the dust grains absorb the stel- 
lar spectrum with respect to the Planck mean opacity Kp(Tt) at a 
temperature r» of the stellar photosphere. In contrast, the FLD 
approximation assumes that dust grains are embedded in a lo- 
cally isotropic radiation field and that their absorption behavior 
is based on the Rosseland mean opacity kj^{T) at the local radia- 
tion temperature T at the top of the cavity. Computing the Planck 
or Rosseland mean opacity for a given temperature T mostly re- 
sults in a difference of only a factor of two. However, computing 
the Planck mean opacity with the stellar temperature instead 
of the local radiation temperature T leads to enormous differ- 
ences. In the next subsection, we determine and illustrate the 
difference between these opacities. 

7.1. Opacities 

For the FLD part of our hybrid radiation transport scheme, we 
compute the Rosseland mean opacities Kg^{T) directly from the 
frequency-dependent opacity table (Laor & Draine 1993). To ob- 
tain the absorption coefficient in a given grid cell, these opacities 
are multiplied by the local dust-to-gas mass ratio Mdu,st/A^gas(J^)- 
In our models, the evaporation temperature of dust grains is cal- 
culated depending on the local gas density by applying the for- 
mula of Isella & Natta (2005), which represents a power-law fit 
to the Pollack et al. (1994) data. For the opacities in the simula- 
tion of Krumholz et al. (2009), the authors use linear regressions 
to the Pollack et al. (1994) data in a small couple of intervals; 
the density dependence of the evaporation of dust grains is not 
taken into account. 

The Rosseland kj^{T) and Planck Kp(T) mean opacities used 
in our simulations and those applied in Krumholz et al. (2007), 
Eq. (11), are visualized in Fig. 6. Owing to the similarities be- 
tween our and the Krumholz et al. opacities, the resulting ra- 
diation pressure in the cavity and shell should be comparable 
in the runs using the FLD radiation transport method. If there 
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were a difference at all, the higher mean opacities in our simu- 
lations for T > 600 K would lead to a slightly higher radiation 
pressure onto the swept-up mass shell than in the simulations of 
Krumholz et al. If our opacities in the FLD run lead to an un- 
stable cavity shell, the opacities used by Krumholz et al. have to 
yield the same outcome. 

In our RT-i-FLD runs, the FLD approximation is only used 
to determine the thermal radiation by dust grains. In addition, 
the irradiation of the stellar luminosity up to the first absorption 
event (i.e. up to an optical depth of t(v) » 1) is computed by a RT 
step that appropriately resembles the propagation and absorption 
of the stellar photons. The dust grains absorb the star light ac- 
cording to the frequency-dependent absorption coefficient k(v). 
Although this frequency dependence is fully covered in our RT 
solver (we use 79 frequency bins in the case of the lOr & Drame 
( i 99 ) opacities), in this derivation - for simplicity - we use 
the gray approximation by computing the corresponding Planck 
mean opacities Kp{Tt) that depend on the stellar surface temper- 
ature r,. Fig. 7 shows the resulting Planck mean opacities Kp{T^,) 
as a function of the stellar surface temperature . To obtain the 
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Fig. 7. Planck mean opacities per gram gas as a function of 
the radiation temperature r» of the stellar photosphere as used in 
the RT step of the hybrid radiation transport scheme in Kuiper 
et al. (2010a, 201 1) and in this paper. The upper horizontal axis 
shows the mass of the proto-star, which corresponds to the stellar 
surface temperature Tt, taken from the losokawa & Omukai 
(200 ) evolutionary track for a constant accretion rate of M = 
10"^ Mo yr-i. 



absorption coefficient Kp(Tt) in a given grid cell of the computa- 
tional domain for these stellar photons, the opacities Kp(Tt) are 
multiplied by the evaporation probability, and, therefore, also de- 
pend on the local gas density and temperature. 

The mass of the proto-star during the onset of the cavity shell 
instability in the FLD run is about 20 Mq and higher. Comparing 
the absorption coefficient of the dust grains in the shell on top of 
the optically thin cavity in the FLD approximation (Fig. 6) with 
the RT method (Fig. 7), immediately illustrates that the FLD ap- 
proximation underestimates the absorption coefficient tremen- 
dously (at least by a factor of 20-50). 

In the following subsections, we show that this underesti- 
mated absorption coefficient for direct stellar irradiation in the 
FLD method leads to epochs of marginal Eddington equilib- 
rium in the cavity shell, allowing the shell to become unsta- 
ble, whereas the higher-order treatment of stellar irradiation via 



RT results in a stable cavity shell with radial velocities up to 
V > 100 km s"' and mass-loss rates of M ^ 2x 10^'^ 



Mo yr 



7.2. Radiative acceleration 



To compute the radiative acceleration arad of the swept-up mate- 
rial at the top of the optically thin cavity, we derive analytic ex- 
pressions for the different radiation transport approaches based 
on formulae similar to the well-known Eddington limit. We start 
with the formula for radiative acceleration 



F 

«rad ^ K — 
C 



(6) 



for the radiative flux F and the speed of light c. Since we only 
wish to estimate the radiation pressure by the first absorption of 
the (isotropic) stellar light and the cavity is treated as optically 
thin up to the location ^cavity of the swept-up material on top of 
it, we compute the radiative acceleration onto dust grains only 
along the direction of the polar axis (ID) and express the radia- 
tive flux by 



^(^cavity) — 



cavity 



(7) 



Eq. 7 is valid either for the assumption of an optically thin cav- 
ity or - owing to the conservation of momentum and energy - 
as long as the swept-up mass in the cavity shell can be treated 
as spherically symmetric. By means of the absorption of stel- 
lar photons, the dust grains heat up and then re-emit photons at 
longer wavelengths. Integrating over all photons penetrating a 
shell at the radius r and the momentum gained by the dust grains 
ensures that there is an excellent conservation of momentum. 

Inserting Eq. 7 into Eq. 6 results in a radiative acceleration 
that depends on the absorption coefficient k, the stellar luminos- 
ity L», and the cavity height ^cavity above the star 



Orad 



AncR^ ., 

cavity 



(8) 



This radiative acceleration has to be compared with the gravita- 
tional attraction Ograv of the dust grains at the cavity shell position 
^cavity of the Central stellar mass M» 



cavity 



(9) 



where G is the gravitational constant. The gravity of the less 
massive circumstellar accretion disk is neglected here. In obser- 
vations and simulations, even the mass of the accretion disk plus 
the rotating large-scale torus further away from the cavity shell 
is less than or comparable to the mass of the central star; hence 
neglecting this mass reduces the stellar gravity by only a factor 
of two at most. To minimize the numbers of independent vari- 
ables, we derive the stellar luminosity L» from the stellar mass 
Mt by assuming that the proto-star follows the evolutionary track 
of ' losokawa & Omukai (20( ) with a constant accretion rate of 
M - 10"^ Mq yr~'. The assumption of this particular evolution- 
ary track does of course influence the exact quantitative compar- 
ison of the radiative versus the gravitational acceleration, but the 
qualitative outcome does not change. The critical factor in this 
comparison is the strength of the absorption coefficient k. 

In the FLD approximation, the absorption coefficient k is 
given by the Rosseland mean opacity kr(T) at the local radia- 
tion temperature T of the cavity shell. 
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In the following, we compute the local gas temperature T 
at the top of the cavity by the formula of Spilzer (1978) for an 
optically thin stellar environment with the prerequisite ^cavity » 

that the distance /^cavity to the shell is much larger than the 
stellar radius 



r(/?cavity) = (0.5 

\ ^cavity / 



(10) 



where jSopadty represents the exponent of the slope of the 
frequency-dependent opacities at longer wavelengths. The as- 
sumption of gray absorption would, e.g., imply ySopacity = and, 
therefore, results in a temperature slope of T oc iJ^T^^j^. The slope 
for the opacities of Laor & Draine (1993) is ySopadty ~ 2.05. 

We can express the radiative acceleration in the FLD case 
using Eq. 8 with the Rosseland mean opacities of Fig. 6 at the 
local radiation temperature T(Rcavity) computed via Eq. 10. The 
radiative acceleration in the RT case is computed by using Eq. 8 
with the corresponding Planck mean opacities (including the lo- 
cal dust evaporation probability) of Fig. 7 at the stellar temper- 
ature Tt, which is for a specific stellar mass M» also taken from 
the evolutionary track of Hosokawa & Omukai (2009). The final 
formulae for the accelerations are given by 



a 



a 



cavity 

L*(M*) 



rad 



a 



grav 



R^ , ■ 

cavity 



(11) 

(12) 
(13) 



The two remaining independent parameters are the actual stellar 
mass M, and the location /?cavity of the cavity shell material. In 
the RT case, the opacity in regions cooler than the dust evapo- 
ration temperature is independent of the location /^cavity, i-C- the 
radiative acceleration scales with a cc R^^ in the same way as 

cavity 

the gravitational acceleration. The ratio af^^/ogmv of these accel- 
erations, therefore, is scale-free and depends only on the proto- 
stellar properties leading to the well-known formula for the gen- 
eralized Eddington limit 



47tGc 



(14) 



7.3. Results 



In the following, we compare the accelerations depending on 
the remaining variables M» and /?cavity Therefore, Fig. 8 com- 
pares the radiative accelerations in both the FLD approxima- 
tion and the RT approach to gravity at two different locations 
from the proto-star as a function of the stellar mass M*. The 
steep increase in the radiative forces at roughly 6-8 Mg marks 
the point in time at which the stellar luminosity in the Hosokawa 
& Omukai (2009) tracks starts to dominate over the accretion 
luminosity. 

The results of the analytical estimate derived in the previ- 
ous section and visualized in Fig. 8 fully support the outcome 
of the radiation hydrodynamical simulations: In the RT case, 
the cavity shell becomes super-Eddington (ratio of radiative to 
gravitational acceleration of unity) for a 10 Mg star and be- 
yond. Afterwards, the radiative acceleration increases rapidly 
to be one or two orders of magnitude higher than gravity. The 



cavity shell is highly super-Eddington in all RTh-FLD simu- 
lations. Furthermore, the continuous enhancement in radiative 
acceleration from a 10 Mg up to a 26 Mg star explains the 
difference between the first and second outflow launch in the 
low-mass (Mcore = 50 Mg) case, which are not detected in the 
^core - 100 Mg cases, yielding higher accretion rates and more 
rapid stellar evolution. 

At large radii, also in the FLD runs, the cavity shell is fi- 
nally super-Eddington, but even for a 30 Mg star only by a fac- 
tor of a few. Most of the time, the cavity shell in the FLD case is 
in marginal Eddington equilibrium. In the runs for the FLD ap- 
proximation, the cavity shell can be estimated to be in marginal 
Eddington equilibrium for a proto-stellar mass of roughly 20 Mg 
and an extent of the cavity of roughly 2000 AU (left panel of 
Fig. 8). As a consequence, the launched cavity shell stops its ex- 



pansion along the polar axis at a height R^.^ 



cavity 



2000 AU, where 



it undergoes an epoch of marginal Eddington equilibrium, i.e. ra- 
diative forces are balanced by gravity. In this situation, these 
cavity shells are supposed to become radiative Rayleigh-Taylor 
unstable (Jacquci ^: Ivrumholz 201 1). This is in full agreement 
with the FLD simulation results herein. 

This behavior changes drastically when the direct absorption 
of stellar light is taken into account. In the simulations, which 
treat the stellar iiTadiation with a RT step, the dust grains of the 
cavity shell absorb the stellar light according to its frequency. 
In this way, the acceleration of the shell turns out to be one to 
two orders of magnitude higher than in the FLD approximation. 
In the analytical estimate of the RT case, the environment of 
the 20 Mg proto-star is super-Eddington by a factor of 20 and 
increases at large radii up to 100 for a proto-star of 30 Mg . 

In the simulations, the difference in the radiation transport 
methods exists only up to an optical depth of t(v) « 1 in the cav- 
ity shell, i.e. up to the location where all stellar photons are ab- 
sorbed. In our hybrid radiation transport scheme, the re-emission 
of the photons by dust grains (mostly at longer wavelengths in 
the infrared regime) is computed in the FLD solver step. The 
corresponding length scale is given by /* - {Kt{v) Pcavity)"' and 
therefore depends on the frequency v of the photons as well as 
on the actual shell density Pcavity This absorption length scale /» 
is shown as a function of frequency v for different densities p of 
the cavity shell in Fig. 2. 



8. Comparisons 

8.1. Comparison with a 3D gray FLD simulation 

Krumholz et al. (2009) presented a self-gravity radiation hydro- 
dynamics simulation of the collapse of a 100 Mg pre-stellar core. 
The outer core radius was chosen to be 0.1 pc and the den- 
sity profile declines in proportion to r The initial isother- 
mal temperature of the core is 20 K. The pre-stellar core is ini- 
tially in solid-body rotation without any turbulent motion. The 
model describes a so-called monolithic collapse scenario. The 
applied radiation transport method is the gray flux-limited Af- 
fusion approximation. These properties of this configuration are 
the same as in our simulation run "C-FLD". The equations are 
solved on a 3D adaptive mesh refinement grid in cartesian coor- 
dinates. Densities above the Jeans density on the finest grid level 
are represented by sink particles. 

During the simulation, bipolar "radiation-filled bubbles" are 
blown into the environment of the forming massive star To clar- 
ify our terminology, these "radiation-filled bubbles" correspond 
to the "radiation-pressure-dominated cavities" and the "bubble 
wall" is called a "cavity shell" in this paper At an extent of 
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roughly from 1200 to 2000 AU (from Krumholz et al. (2009), 
Fig. 1 therein) the cavity shell is subject to a "radiative Rayleigh- 
Taylor instability", meaning that the radiation pressure expands 
the optically thin region only at specific solid angles, while at 
other solid angles the concentrated mass load on top of the cav- 
ity shell is able to penetrate into the low-density region. 

This mechanism of shell instability, the focus of radiation 
pressure onto solid angles with lower optical depth, and the col- 
lapse of condensed material from the infalling envelope at other 
solid angles, precisely matches the outcome of our simulations, 
if the FLD approximation is used for the stellar radiation feed- 
back. After the epoch of instability - the caving-in of material 
- the continuing evolution of their simulation differs in terms of 
morphology from our simulations because of the evolving non- 
axisymmetric structures. We comment further on this difference 
after the following paragraph on the RT + FLD simulations. 

In our simulations using the frequency-dependent RT step 
for the stellar irradiation, the mass load on top of the cavity 
shell causes a dependence of the optical depth on the solid an- 
gle (owing to the decreasing centrifugal forces towards the pole). 
However, in these simulations, the caving-in of material at solid 
angles of high optical depth is prohibited; this is most likely 
caused by the treatment of the direct stellar irradiation by means 
of a RT approach preserving the natural isotropy of the irradi- 
ation up to the first absorption layer. Furthermore, the RT ap- 
proach accounts for the respective frequency-dependent absorp- 
tion coefficients as demonstrated in the previous sections. 

The cavity is affected by Rayleigh-Taylor-like instabilities, 
when the FLD approximation is used for the radiation transport. 
This implication is verified by the axial symmetric simulations 
performed herein as well as by the 3D simulation in Krumholz 
et al. (2009). Improving the radiation transport scheme by in- 
corporating a ray-tracing of the direct stellar irradiation corrects 
for this behavior and leads to a stable cavity growth. This impli- 
cation is verified by the axial symmetric simulations performed 
herein as well as by the 3D simulation in Kuiper et al. (20! I). 
The largest differences between the numerical treatment in the 
simulation of Krumholz et al. (2009) and the simulations with 
similar initial conditions presented herein are twofold: the radi- 
ation transport approaches and our assumed axial symmetry. As 
mentioned above, the simplification to axial symmetry increases 
the difficulty of a comparison of epochs after the occurrence of 
the instability, but the overall simulation results indicate that the 
outcome, regardless of whether the instability occurs, does not 
depend on the axial symmetric approach: 

- The 2D simulations in the FLD approximation match the 
outcome of the 3D simulation in the FLD approximation, 
namely the instability of the cavity. 

- In the 2D simulations with the improved radiation transport 
scheme (ray-tracing -i- FLD), the cavity remains stable. 

- In our 3D simulation (Kuiper el al. 201 using a frequency- 
dependent RT step to account for stellar irradiation, the out- 
flow cavity contains strong non-axisymmetric features, but 
no instability is detected. 

This implies that the stability of the outflow cavity is strongly 
effected by the direct stellar irradiation of the massive star, inde- 
pendent of the symmetry. 

8.2. Comparison with frequency-dependent FLD simulations 

In Yorke & Sonnhalter (2002), the authors presented six self- 
gravity radiation hydrodynamics simulations of the collapse of 
a pre-stellar core. The initial core mass in the simulations was 



chosen to be 30, 60, and 120 Mq with an outer core radius of 
0.05, 0.1, and 0.2 pc. The density profile drops proportional to 
r^^. Since in these simulations an additional inflow of mass at 
the outer boundary into the computational domain is allowed, the 
total mass reservoir of the forming star is not limited to these ini- 
tial core masses. The initial isothermal temperature of the core is 
20 K. The pre-stellar core is initially in solid-body rotation with- 
out any turbulent motion. The equations are solved on a static 
nested grid of three levels in cylindrical coordinates assuming 
axial symmetry and midplane symmetry. The forming massive 
star is represented by a sink cell at the origin of the domain. 
The radiation transport method used is the gray as well as the 
frequency-dependent FLD approximation. 

The outflow properties in the simulations with gray FLD and 
initial core masses of 30 and 120 Mq are neither visualized nor 
discussed in Yorke & Sonnhalter (2002). In the simulation with 
gray FLD and an initial core mass of 60 M© , the massive star 
stops its mass growth at = 20.7 M© and maintains a lumi- 
nosity of Lt = 5.2 X 10"* Lq. No polar cavity forms before the 
end of the simulation at 1 10 kyr. 

In the three simulations with frequency-dependent FLD, po- 
lar cavities are formed. In the 30 Mq case, the cavity interacts 
with the infalling envelope and collapses within a few thousand 
years. In the two higher mass cases, the infall even at larger 
radii is reversed shortly after the launch of the outflow cav- 
ity; therefore the cavity does not interact with a gravitationally 
driven infall anymore and simply expands in time. The 60 Mq 
simulation run was stopped after an evolutionary age of 45 kyr. 
Nevertheless, in the 120 Mq case, the radiation pressure eventu- 
ally evacuates the stellar environment in all directions, including 
the disk midplane. After its first expansion phase, this radiation- 
pressure-dominated cocoon-like structure decreases again, espe- 
cially in the direction of the highest optical depth (the midplane), 
similar to the behavior of the expansion of the cavity shell in 
the FLD simulation by Krumholz el al. (200'-') as well as in 
our FLD simulations. The expansion velocity of the cavities of 
V > lOkms ' corresponds more to the velocities observed in 
our FLD simulation (v = 6 - 30 km s"') than those using the 
ray-tracing method (v ~ 100 km s"'). 

The frequency dependence of the FLD solver (as in our RT 
method) might also be important. Using the FLD approximation 
in the gray limit leads to a downgrading of the radiation field 
throughout the cavity region, i.e. the dust grains in the cavity 
shell absorb the stellar flux with the Rosseland mean opacity at 
the local radiation temperature of the cavity shell, which does 
not resemble the broad stellar irradiation spectrum. By dividing 
the stellar flux into several frequency bins and treating the flux of 
each bin in the FLD approximation, one also can ensure that the 
high-frequency part of the stellar spectrum is absorbed with the 
appropriate absorption coefficient. In addition to a more correct 
description of the absorption behavior in the cavity shell, the RT 
approach provides a most realistic representation of the photon 
paths emitted at the stellar photosphere. Owing to the integral 
over the solid angle in the FLD approximation, the resulting flux 
does not include the correct propagation direction. 

Moreover, including the long-wavelength part of the stellar 
spectrum explicitly can accelerate the layers on top of the swept- 
up cavity shell, which is optically thick for the short wavelengths 
only. Furthermore, Schartmann et al. (201 1) demonstrated that 
in numerical simulations using a RT approach for the radiation 
feedback, the acceleration of a gas clump by radiation pressure 
against a gravitational potential is inversely proportional to the 
optical depth of the clump, as expected by analytic estimates 
similar to the computation of the Eddington limit. 
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8.3. Comparison witli an analytic stability analysis 

Jacquet & Kmmholz (2011) derived a criterion for the stabil- 
ity of radiation-pressure-dominated cavities for the adiabatic 
approximation. They provided growth times of the radiative 
Rayleigh-Taylor instability in cavity shells around massive stars. 

This analytic analysis fully supports the outcome of our in- 
vestigations. From our estimate of the accelerations in Sect. 7, 
we compute the expansion timescale of the cavity shell to be 

J 2 ^cavity /i c\ 

■ (Ij) 
^rad ^grav 

In the case of FLD simulations, the radiation pressure onto the 
cavity shell is in marginal equilibrium with gravity {E ^ 1) 
for long epochs in time. The corresponding shell expansion 
timescale goes to infinity during these epochs, allowing the ra- 
diative Rayleigh-Taylor instability to set in. In the case of RT 
simulations, radiation pressure exceeds gravity by 1-2 orders 
of magnitude, leading to a much shorter cavity shell expan- 
sion timescale. For a stellar mass of 20 M© and a cavity radius 
of ^cavity = 1400 AU, One obtains an expansion timescale of 
fcavity ~ 0.68 kyr. 

Jacquet & Krumholz (201 1) stated that the growth times of 
the radiative Rayleigh-Taylor instability in cavity shells around 
massive stars should always be shorter then 1 kyr, but actually 
the numbers one reads from their plots (Fig. 3 and 4) are in- 
stead 6.3 - 40 kyr depending on the actual stellar mass and the 
wavelength of the instability. Hence, the expansion timescale 
of the shell is at least one order of magnitude shorter than the 
timescale for the growth of the radiative Rayleigh-Taylor insta- 
bility. Even, for the largest shell radius of /^cavity - 0.1 pc at 
the outer boundary of our computational domain, the expansion 
timescale is fcavity ~ 10 kyr During this 10 kyr, the massive star 
increases its mass approximately from 20 to 30 by disk ac- 
cretion. Therefore, a cavity shell that has formed around a mas- 
sive star is expected to be able to expand up to the stellar cluster 
scale (^ 0.1 pc) before being prone to the radiative Rayleigh- 
Taylor instability. 

One might discuss whether an actual outflow of a massive 
star might fulfill the conditions for adiabaticity. In Jacquet & 
Krumholz (20 1 i), the authors assume that the cavity shell is in 
the "adiabatic regime". This is estimated in Jacquet & Krumholz 
(2011) in Eq. (84) under the assumption of a hot cavity shell 
with T < 1100 K at a very large radius of ^cavity ~ 10"* AU. 
In the simulation data of Krumliolz et al. (2009) and our simu- 
lations in this paper, we observe much lower shell temperatures, 
even at much smaller radii. As a consequence, for the stellar evo- 
lutionary tracks of Hosokawa & Omukai (200 ) and the shell 
temperature estimate given in Eq. (10), the cavity shell is in the 
non-adiabatic regime even out to an extent of 0. 1 pc for a stellar 
mass of Mt < 40 M©. In the non-adiabatic regime, the radiative 
Rayleigh-Taylor instability is damped by diffusion. The quanti- 
tative change caused by the damping remains unclear, and the 
growth rates of the radiative Rayleigh-Taylor instability in the 
non-adiabatic regime have not yet been derived. 

8.4. Comparison with observations 

As mentioned in the introduction, radiation-pressure-dominated 
cleared cavities or optically thin bubbles of the discussed sizes 
around massive proto-stars have not yet been observed. Similar 
cavities are also proposed to form during a non-radiative but 



magnetized pre-stellar core collapse as shown in MHD simu- 
lations of Banerjee & Pudrilz (2006, 2007) and Hennebelie el al. 
( ). The observation of cavities, especially in their initial 
launch phase, are hampered by the extinction (the cavities are 
still embedded in the infalling envelope) as well as the lim- 
ited resolution of the inner core structure close to the proto-star. 
However, large-scale polar cavities could potentially diminish 
the extinction of the envelope at pole-on viewing angles. 

Moreover, a direct comparison of the detailed morphology 
of the simulated outflow cavities is of course restricted by the 
limitations and caveats of the numerical simulations discussed 
in Sect. 2. Nevertheless, the basic properties of the radiation- 
pressure-dominated outflows investigated in this paper include 
a shell velocity on the order of v a; 100 km s as well as a 
mass-loss rate of roughly M x IQ^'^ Mq yr"' (see Fig. 5). 

9. Summary 

We have performed six self-gravity radiation hydrodynamics 
simulations of a collapsing pre-stellar core including the launch 
of a radiation-pressure-dominated outflow cavity by the centrally 
forming star For the radiation transport method, three of these 
simulations were performed using the FLD approximation and 
three use the hybrid radiation transport scheme introduced in 
Kuiper et al. (20 It) ). In the hybrid radiation transport scheme, 
the stellar irradiation is computed in a frequency-dependent RT 
step and only the computation of the thermal dust emission 
makes use of the FLD approximation. The three different ini- 
tial conditions of the simulations vary in terms of the initial pre- 
stellar core mass (50 and 100 Mq) as well as the initial density 
slope (p oc r"'-^ and p cc r^^). The analysis of these simula- 
tions indicates that only in the case of the FLD approximation 
does the cavity shell undergo an epoch of marginal Eddington 
equilibrium, i.e. radiative forces are balanced by gravity. As a 
consequence, the expansion of the outflow cavity along the po- 
lar axis stops and the infalling envelope material gathered on top 
of the cavity shell caves in. In contrast to this scenario, simu- 
lations using the hybrid radiation transport scheme predict that 
there are stable and rapidly growing outflow cavities, and that no 
long-term epochs of marginal Eddington equilibrium occur. 

We analyzed the outcome of these simulations in more de- 
tail using analytical estimates of the radiative acceleration and 
the growth timescale of the cavities with respect to the radiation 
transport method. In contrast to the FLD assumption that dust 
grains absorb photons according to the local radiation tempera- 
ture, dust grains in the shell on top of the optically thin cavity are 
directly irradiated by photons from the luminous massive star. 
Therefore, including the feedback effect of direct stellar irradi- 
ation via a frequency-dependent RT approach leads to the pre- 
diction of stronger radiative forces in the cavity shell. Including 
the direct absorption of stellar photons enhances the accelera- 
tion of the shell by between one and two orders of magnitude 
in comparison with the FLD assumption. This analytic estimate 
confirms the much higher shell velocity detected in the radia- 
tion hydrodynamics simulations including the RT step for stellar 
radiation feedback. These estimates of the relevant timescales 
imply that this enhancement in acceleration leads to a cavity ex- 
pansion time that is much shorter than the timescale needed for 
the radiative Rayleigh-Taylor instability to occur 

The results of this study does not prohibit in general the exis- 
tence of (radiative) instabilities in cavity shells. The simulation 
results for the three different initial conditions show a depen- 
dence of the radiative launching and cavity growth on the den- 
sity distribution of the stellar environment, e.g. initially flatter 
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density profiles lead to larger amounts of mass at larger radii 
and, therefore, a heavier mass load on top of the cavity shell. 
The limitations and caveats of the numerical simulations (e.g. the 
assumptions of no magnetic fields and perfect dust-to-gas cou- 
pling) do not allow a complete or even comprehensive descrip- 
tion of the complex outflow region of a massive star. 

Nevertheless, this study clearly emphasizes the strong influ- 
ence of the (spectral) stellar irradiation on both the dynamics 
and the stability of the first absorption layer. It has been demon- 
strated that a careful numerical treatment of this direct stellar 
irradiation feedback can be achieved by using a RT approach. 
For comparison with future high-resolution observations of the 
inner core regions, the velocity and the mass-loss rate of these 
radiation-pressure-dominated outflows are given: in the simula- 
tions, we find a shell velocity on the order of v » 100 km s"^ and 
a mass-loss rate of roughly M » 10""' M© yr"^ The qualitative 
and quantitative analyses of our numerical simulations, includ- 
ing different initial conditions and radiation transport methods 
combined with the analytical estimates of the relevant forces and 
timescales, cast severe doubt on the occurrence of the radiative 
Rayleigh-Taylor instabiUty in the radiation-pressure-dominated 
cavities around forming high-mass stars. 
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Fig. 3. Simulation snapshots of the gas density for one of the FLD (left panels) and one of the ray-tracing + FLD (right panels) runs 
for three different points in time during the launch of the radiation pressure dominated cavities. The data is taken from the runs with 
initial core mass Mcore = 50 and a radial density slope of /3 - -2. 

Note: The spatial section of the RT+FLD snapshots increases with time to follow the rapid expansion of the outflow cavity. 
Animations of the launch of the radiation-pressure-dominated cavities in the simulations are available as online material, too. 
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Fig. 4. Size of the cavity /?cavity as a function of time t for the three different initial conditions. All runs were performed using the 
RT-hFLD hybrid radiation transport scheme (solid lines) and the FLD scheme (dashed lines). The lowest and highest extent of /?cavity 
of 10 AU and 0.1 pc, respectively, are given by the inner and outer rim of the computational domain. 
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Fig. 5. Gas density p (upper left panel), temperature T (upper right panel), radial velocity v (lower left panel), and radial mass-loss 
rate (lower right panel) as a function of height z above the star along the polar axis. Data is for the run with initial core mass 
McoK = 100 Mq and radial density slope /? = -2. The snapshot in time (f = 16.7 kyr) is chosen in such a way that the cavity shell 



has arrived at the same location in both simulations (RT-hFLD and FLD only). The lowest density of p = 10 
left panel is given by the numerical floor value of the hydrodynamics solver. 
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Fig. 6. Rosseland mean opacities atr (left panel) and Planck mean opacities Kp (right panel) per gram gas as a function of the local 
radiation temperature T. The label "FLD 1" denotes the opacity description by Krumholz et al. (2007), Eq. (11) therein. The label 
"FLD 2" denotes the opacities used in Kuiper et al. (2010a, 201 1) as well as in this paper; these opacities are computed based on 
the density-dependent evaporation of dust grains using the low density p = 10 g cm"^ of the polar cavities. 
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Fig. 8. Radiative and gravitational acceleration at a position 7?cavity = 2000 AU (left panel) and Rcavky - 10000 AU (right panel) 
above the massive proto-star as a function of the stellar mass M, for three different radiation transport approaches: The label "FLD 
1" denotes the FLD approximation with the opacity description of Krumliolz el al. (2007). The label "FLD 2" refers to the FLD 
approximation with the opacities used in Kuiper et al. (2010a, 201 1) and in this paper. The label "Ray-Tracing" refers to the RT step 
as in the hybrid radiation transport scheme. 



